%% ARP efficiency calibration 2023/01/04
t = [0 0.1000 0.2500 0.5000 1 2 4 5]; % [ms]
pkOD = [112.1161 210.6172 338.9238 501.1355 708.7992 930.6653 1.0355e+03 1.0727e+03];
pkOD_err = [11.7481 10.2786 10.4597 20.3377 17.9409 49.5823 36.4666 44.7537]; % standard deviation of five measurements

% subtract background (counts at t=0)
pkOD = pkOD - pkOD(1);
pkOD_err = sqrt(pkOD_err.^2 + pkOD_err(1).^2);

% get the atom number ratio with errorbar. Assume at t=4, atom number is 1
y = pkOD; 
e_y = pkOD_err;
yref = pkOD(7);
e_ref = pkOD_err(7);

ratio = y./yref;
err_ratio = ratio.*sqrt((e_y./y).^2+(e_ref./yref).^2);

%% Number dependence
h=figure(172)
ratio = [0 0.1067 0.2456 0.4213 0.6462 0.8865 1.0000];
err_ratio = [0 0.0175 0.0199 0.0309 0.0355 0.0663 0.0587];

% N=2, Q2
% data: '20230416_100VR_30G_10kHz_REMPI_scan_7' 
% Cyclescut = {[]}
y2 = [1.0000    7.1000   16.2000   23.7000   43.8000   56.8000   63.8000];
y2 = y2-y2(1);
err_y2 = [0.5477    0.8888    1.3115    1.5906    2.1260    2.3958    2.5338];


% N=1 
% data:'20230416_100VR_30G_10kHz_REMPI_scan_8'
% cyclesCut =  {[]}
y1 = [2.5000    2.7500    7.4167   12.2500   21.0000   23.0833   26.5833];
y1 = y1-y1(1);
err_y1 = [0.5270    0.5590    0.8620    1.0308    1.3281    1.4410    1.5069];

% N=0
% data: '20230416_100VR_30G_10kHz_REMPI_scan_9'
% cyclesCut = {[]}
y0 = [4.9000    9.1000   10.8000   15.4000   19.6000   24.9000   28.7000];
y0 = y0-y0(1);
err_y0 = [0.7416    1.0149    1.0677    1.2570    1.4422    1.6093    1.7117];

colors={[0 0.4470 0.7410],[0.8500 0.3250 0.0980],[0.9290 0.6940 0.1250],[0.4940 0.1840 0.5560],[0.4660 0.6740 0.1880],[0.3010 0.7450 0.9330],[0.6350 0.0780 0.1840]};

errorbar(ratio,y2,err_y2,err_y2,err_ratio,err_ratio,'o','markersize',8,'linewidth',2,'displayname','N=2','color',colors{1});
hold on
errorbar(ratio,y1,err_y1,err_y1,err_ratio,err_ratio,'^','markersize',8,'linewidth',2,'displayname','N=1','color',colors{2});
errorbar(ratio,y0,err_y0,err_y0,err_ratio,err_ratio,'s','markersize',8,'linewidth',2,'displayname','N=0','color',colors{3});

% Linear fit
ft = fittype( 'poly1' );
[fitresult2, gof] = fit(ratio', y2', ft );
[fitresult1, gof] = fit(ratio', y1', ft );
[fitresult0, gof] = fit(ratio', y0', ft );
xx = linspace(0,1.05,1000);
yy2 = fitresult2(xx);
yy1 = fitresult1(xx);
yy0 = fitresult0(xx);

plot(xx,yy2,'-','linewidth',1.5,'color',colors{1})
plot(xx,yy1,'-','linewidth',1.5,'color',colors{2})
plot(xx,yy0,'-','linewidth',1.5,'color',colors{3})

xlabel('Normalized atom number');
ylabel('KRb ion counts')
box on; %grid on; %grid minor;
set(gca,'Fontsize',16,'linewidth',1.5)

legend('N=2','N=1','N=0')
h.Position = [100 100 600 400];
xlim([0 1.08])
ylim([-5 70])
% extract branching ratio from fitted slope
y = [53.86 38.1 24.81];
e = [(62.01-53.86)/2 (45.14-38.1)/2 (30.58-24.81)/2];

tot = sum(y);
e_tot = sum(sqrt(e.^2));
ratio = y/tot;
e_ratio = y/tot.*sqrt((e./y).^2+(e_tot/tot).^2);

%% Fit ARP to Landau-Zener formula
[xData, yData] = prepareCurveData( t, pkOD );

% Set up fittype and options.
ft = fittype( 'a*exp(-b*x)+c', 'independent', 'x', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Lower = [-Inf 0.1 -Inf];
opts.StartPoint = [0.99242382885004 0.432017372985733 0.501694876211973];

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );

% Plot fit with data.
figure( 'Name', 'Exponential fit' );
h = plot( fitresult, xData, yData);
legend( h, 'pkOD vs. t', 'Exponential fit', 'Location', 'NorthEast' );
set(gca,'fontsize',16)
set(gcf,'color','white')

% Label axes
xlabel t
ylabel pkOD
grid on